source('00_util_scripts/mod_bulk.R')
library(clusterProfiler)

tcc7a <- read_delim('mission/Duan_GSDMx/ttc7a.mut.bulk.csv')

tcc7a |>
  filter(gene == 'TMEM33') |>
  pivot_longer(BCH_disease:VPS6) |>
  mutate(group = ifelse(name == 'BCH_disease', name, 'Ctrl')) |>
  ggplot(aes(group, value)) +
  geom_jitter(height = 0, width = .01) +
  stat_summary(geom = 'crossbar', fun = mean, width = .3) +
  annotate('text', x = 1.5, y = 12, label = 'adjusted P = 0.343') +
  labs(title = 'TMEM33 expression', y = 'Normalized expression') +
  theme_pubr() +
  expand_limits(y = 0)

# GO gsea ---------
tcc7a |>
  filter(pvalue < .05) |>
  count(log2FoldChange > 2)

tcc7a.gsego <- tcc7a |>
  filter(padj < .05, abs(log2FoldChange) > 2) |>
  pull(log2FoldChange, name = id) |>
  sort(decreasing = T) |>
  gseGO(ont = 'BP',OrgDb = 'org.Hs.eg.db',keyType = 'ENSEMBL',pvalueCutoff = 1) 

tcc7a.gsego@result |>
  as_tibble() |>
  filter(str_detect(Description, 'culum')) |>
  select(Description, NES, p.adjust) |>
  slice_max(abs(NES), n = 10) |>
  mutate(Description = str_wrap(Description, 40) |> fct_reorder(NES)) |>
  ggplot(aes(NES, Description, fill = p.adjust)) +
  geom_col() +
  theme_pubr(legend = 'right') +
  scale_fill_gradient(low = 'red', high = 'blue')

tcc7a.gsego |>
  dotplot()

tcc7a.gsego@result |>
  filter(str_detect(Description, 'phagocytosis'))

go.bca <- read_delim('mission/Duan_GSDMx/msigdb.manual.txt', delim = ';') |>
  separate_longer_delim(gene, ',') |>
  pull(gene)

tcc7a |>
  filter(gene %in% go.bca, padj < .05) |>
  write_csv('B.cell.activate.BCH.downregulated.GO0042113.csv')

library(enrichplot)
gseaplot2(tcc7a.gsego, 'GO:0034976',
          title = 'GO pathway: response to endoplasmic reticulum stress',
          pvalue_table = T)

# (not) KEGG -------
reactome.iri <- read_csv('mission/Duan_GSDMx/reactome.R-HSA-198933.txt') |>
  separate(MoleculeName, into = c('uniprot','gene'), sep = ' ') |>
  pull(gene)

tcc7a.gseke <- tcc7a |>
  filter(padj < .05, gene %in% reactome.iri)

tcc7a.gseke |>
  ggplot(aes(rank(log2FoldChange), log2FoldChange)) +
  geom_point()

tcc7a.gseke |>
  select(-c(id, baseMean, lfcSE, stat)) |>
  relocate(gene, log2FoldChange, pvalue, padj) |>
  write_csv('Immunoregulatory.BCH.downregulated.R-HSA-198933.csv')
